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"Matched11  Polynomial  Least-Squares  Fitting 
and 

Application  to  Real-Time  Ballistic  Coefficient  Estimation 


Summary 

3 

A  form  of  least-squares  "spline -function"  fitting  is  presented  where, 
in  essence,  a  distinctly  non  -polynomial  function  defined  over  some  gross  inter¬ 
val  is  approximated  by  single  polynomial  expressions  of  degree  M  over  each 
of  contiguous  sub-intervals  with  continuity  conditions  on  the  function  up  through 
the  (M-l)1  st  derivative  being  imposed  at  the  junction  points  (Fig.  1).  Weighted 
least-squares  fitting  to  given  observational  (radar)  data  then  leads  to  linear 
algebraic  equations  for  the  parameter  estimates  instead  of  non-linear  equations 
requiring  solution  by  iteration.  This  technique  leads  to  Eqs.  (4),  (8),  and  (9) 
of  the  analysis  section.  These  relationships  can  be  used  for  general  weighted 
least-squares  fitting  of  observational  data  where  the  "ordinary"  polynomial 
approximation  is  not  sufficiently  accurate.  Application  of  these  relationships 
to  real-time  (3  estimation  appears  promising  and  a  procedure  utilizing  them 
is  outlined  in  this  paper. 

Introduction 

Standard  lea  st  -  squares  polynomial  fitting  methods  are  easy  to  implement 
because  the  parameter  estimates  then  satisfy  linear  algebraic  equations  which 
can  be  readily  solved  by  standard  matrix  inversion  techniques.  However,  long 
segments  of  re-entry  trajectory  are  not  generally  well  approximated  by  a  single 
quadratic  (or  higher  degree)  polynomial  expression.  This  is  basically  the  reason 
why  numerical  integration  of  the  differential  equations  of  the  trajectory  followed 
by  an  iteration  procedure  and  the  ensuing  mathematical  complications  inherent 
in  a  "Maximum-likelihood"  approach,  is  frequently  required  in  order  to  obtain  high- 


*  Ref.  1  discusses  an  alternate  real-time  method. 
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precision  answers.  With  a  view  toward  real-time  applications,  we  there¬ 
fore  attempt  to  circumvent  these  complications  by  finding  a  more  accurate 
functional  representation  than  a  single  polynomial  over  the  gross  interval. 

This  is  accomplished  by  subdividing  the  entire  interval  shown  in 
Fig.  1  into  equal  sections,  assuming  a  different  polynomial  expression  of 
degree  M  over  each  section,  and  requiring  the  continuity  of  the  function  up 
through  the  (M-l)f  st  derivative  at  the  intersection  points  T^,  T  y  ...»  Tj, 

.  .  .  ,  Tp  .  This  nmatchedn  polynomial  approach  will  yield  accuracies  some¬ 
where  intermediate  between  ordinary  least-squares  fitting  and  the  optimum 
iterative  maximum-likelihood  technique,  and  with  considerably  less  compu¬ 
tation  time.  We  now  derive  the  necessary  relationships  and  outline  a  real- 
❖ 

time  technique  that  utilizes  them. 


Analysis: 

In  any  region  --  the  J’th  say  --of  Fig.  1,  a  continuous  differentiable 
function  1  !S(t),f  can  be  approximated  by  an  Mth  degree  polynomial  viz 

M  S<m> 

S<V=  En^r<trTj>m  £orTjstiSTj+i  (1> 

m=0  m  ! 


J  =  1,  2,  .  .  .  ,  p 

where  S.  is  the  value  of  the  function  S(t)  at  t  =  T  in  that  region  and 
Sj^  >Sj^  »  Sj^  .  .  .  are  the  values  of  S(t),  S(t),  and  S(t)  .  .  .  respectively 
at  t  =  Tj  . 


The  "optimum"  real-time  technique  we  define  to  be  "that  technique 
which  provides  parameter  estimates  of  1  SUFFICIENT*  accuracy  for  the 
applications  intended  with  a  minimum  of  computation  time  and  storage". 
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Fig,  1,  Subregions  for  which  different  M'th-degree  polynomial 
expressions  apply. 


A  different  M'th-degree  polynomial  of  the  type  (1)  applies  over  any  region 
other  than  J  .  The  data  times  t^  need  not  be  equally  spaced  in  what  follows 
however,  it  will  be  assumed  that  each  region  has  an  equal  time  span  so  that 
h  =  Tj  -  Tj  ^  is  independent  of  J  . 

Continuity  of  S(t)  up  through  the  (M-l)'  st  derivative  at  the  inter¬ 
section  point  Tj  of  region  J  and  (J-l)  implies  that 


m=p 


(m)  m-|i 
bJ-l  h 


(m-n)! 


U  =  0,  1, 


(M-l) 


(2) 


The  recursion  relationship  (2)  can  be  shown  to  give 


(M-^)  (M-*) 

bJ  "bl 


J-l 

+  E 

k=  1 


v  ,  m 
E  h 


m  =  1 


m: 


c  (M+m-y) 
bk 


v  =  1,  2, 

J  =  2,  3, 


M  (3a) 

P 


Without  the  constraining  relationships  (3),  there  are  (M+l)  p  independent 
parameters  to  be  determined  in  Eq.  (1)  by  the  least-squares  fitting  of  noisy 
S(t)  data.  However,  the  constraints  (3)  show  that  only 

S2(m> . S  <M)  are  independent  parameters  (M  +  p  of  them)  to  be  de¬ 

termined  by  the  fitting.  s/°),  S  (  ^ S  (M-D  are  the  initial  conditions 
at  t  =  Tj  (corresponding  to  the  initiation  of  a  re-entry  track)  and  the 
K  =  1,  2,  .  .  .  ,  p  are  the  M'th  derivatives  of  the  track  in  each  region. 

.It  can  be  shown  by  letting  v  successively  equal  1,  2,  .  .  .  ,  that  Eq.  (3a) 
finally  yields 


(M-v) 

bJ 


v-\ 

=  E  ■ 

k=0 


C(J-l)h] 

k! 


k 


q(M  -v  +k) 
bl 


J-l 

E 

k=l 


[(J-k)1'  -(J-k-l)l']^/t) 


=  1,  2,  .  .  .  ,  M 
J  =  2,  3,  .  .  . ,  p 


(3b) 


In  obtaining  Eq.  (3b),  the  following  relationship  was  utilized: 


J-l  J-l 

E  E 

k=l  k  =1 


J-l 

E 

k  =1 
n 


L'k  <k  ,6lc  <k 

n  n-1  n-1  n-2 


Jkx  <k 


]Sk(M) 

n 


j-i 

E  ~  (J-k-1)  (J-k-2) 


(J-k-n)  4^ 
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where,  for  example, 


k1  <  k 


1  for  k^  <  k 


0  otherwise 


Equations  (3b)  give  the  explicit  dependence  of  the  dependent  parameters 
c  ( M  -  v ) 

,  J  =  2,  3,  .  .  .  ,  p  ;  v  =  1,  2,  .  .  .  ,  M  upon  the  independent  parameters 
The  p'  th  derivative  of  Eq.  (1)  at  t  =  T  can  now  be  written 


1 


M-l  -q 


k=0 


(k+q) 


[MJ-l+e.j)]  hM-\x  J-l 


2  S 


(M) 


(J-k  +  0. 


iJ 


1) 


M-u 


(M-q)!  k=1  k 
<h6iJ>M'" 


{(J- 


(M-q)! 


,(M) 


k  +  6iJ> 


(4) 


M-q 


t.  -  T, 

eiJ s  (-v1) 


Tt  ^  t.^  Tt 
J  l  J  +  l 

J  =  1,  2,  .  .  .  ,  p 

q  =  0,  1,  .  .  .  ,  M 


where  the  following  relationship  was  noted: 


M-l 

2 

k=q 


k  M-k 

x  y 


k!  (M-k)! 


nr  -TT  {<*+y>' 


Notice  that  the  only  restrictions  placed  on  the  arbitrary  function  S(t) 
were  continuity  and  differentiability.  Equation  (4)  is  the  desired  analytical 
expression  for  the  \a' th  derivative  where  |j  *  0  corresponds  to  S(t)  itself. 
As  is  seen  from  Eq.  (4),  the  function  is  represented  over  each  of  p  con¬ 
tiguous  segments  by  a  different  polynomial  of  the  same  degree  M  and  con¬ 
tinuity  conditions  on  S(t)  up  through  the  (M  -  l)1  st  derivative  are  satisfied 
at  all  of  the  junction  points. 

From  Eq.  (4)  we  observe  that 


M-l 

S(t  )  =  2  S 

k=0 


(k) 


[htf-i+e.j)] 


UM  J-l  /  \  >r  \ 

_  .  2l_  s  c  (m) 

k!  +  M!  Sk 

k=  1 


{(j 


„  M 

■k+6ij)  “ 


(J-k+0ij-1)M}  + 


<heij> 

M! 


M 


,(M) 


(5a) 
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as(t.)  [MJ-i+e.j)] 


M* 


1  _ 


(J  =  0,  1,  2,  .  .  .  ,  (M-l)  (5b) 


SS(t.) 

x\ 


MI 


^K,  J8iJ  +  6K<J 


[(J- 


K  +  9iJ) 


M 


(J-K-l  +0  ) 


M 


(5c) 


J  >  K  1  j  2 ,  •  •  •  j  p 


Tt  *  t.  *  T_xl 
J  i  J  +  l 


r  ❖  5jc 

Consider  the  quantity  lS(t.)  -  S.,  J  where  is  given  noisy  experimental 

data  at  t  =  t^  with  variance  c r?  .  Physically,  this  quantity  represents  the 
deviation  of  the  data  from  the  calculated  function  at  time  t^  in  region  J  .  We 
therefore  wish  to  minimize 

N. 


P  J 
A  =  E  E 

J  =  1  i=Nj_  X  +  1 


5jc  2 

S(t.)  -  S. 

V  Y  1 

cr. 
i 


(6) 


where  Nj  is  the  total  number  of  points  up  to  the  value  +  j  . 

Differentiating  A  partially  with  respect  to  S  ^  , 

S^\  K  =  1,  2,  .  .  .  ,  p  and  equating  the  results  to  zero  gives 


Nt 
P  J 

E  E 

J=1  i=Nj  _  1  +1 


S(t.)  -S. 
x  r  l 


cr. 

l 


M-! 


=  0  (7a) 

H  =  0,  1,  ....  (M  -  1) 


P 

2 


nj 

2 


J=1  i=N  +  1 

J  “  1 


S(t.)  -  s 
1  1 


(J. 

1 


* 


8  eM  + 

M!  K,  J  iJ 


+  SK<j[(J-K+VM-(J-K+eiJ-1)M]|=° 

"f  Where  6^7  T  is  the  Kronecker  delta. 

.K,  J 


K  =  1 ,  .  .  .  ,  p 


(7b) 
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Equations  (7)  represent  (M  +  p)  simultaneous  linear  algebraic  equations 
for  the  (M  +  p)  independent  parameter  values.  In  matrix  notation,  Eqs.  (7) 
become 

AX  =  B  (8a) 


where  A  turns  out  to  be  symmetric  and 


A  = 


11 

a  1 2 . 

1 ,  p+M 

21 

a22 . 

"  a2,  p+M 

ap+M,  1  ap+M,  2 . a 


p+M,  p+M 


s 

bl 

l 

bi 

e  (1) 

bl  . 

1 

b2 

X  = 

1 

(M-  1) 

b] 

s 

1  . 

« 

B  = 

1 

1 

1 

1 

I 

1 

s  (M) 

P 

1 

bp+M 

—  m 

(8b) 


There  remains  only  the  determination  of  a  ^  and  b  q.  8=1.  2,  ...  . 

J  a,  p  a 

(p+M)  . 

Equations  (7)  and  (5a)  give 


U+K  .rs- 

J  1  0" 

i 


H,  K  =  0,  1  .  .  .  ,  M- 1  (9a) 


and  aT.  ,  ,  =  a  .  ■ 

K+ 1 ,  1  |jl+  1 ,  K+ 1 


M+pi 
h _ 

V+l,  M+K  "  M! 


2  2 
J  i 


<J-1+V 


M- 


c r. 

l 


SK<J  | 


(9b) 


\i  =  0,  1 ,  .  .  .  ,  M-  1 
K  —1,  2,  p 


and  aM+K,  n+ 1  auL+ 1 ,  M+K* 
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M+K,  M+L 


eM+  a 

(M!)2  II  1  K.J  iJ  K<J 


[<J-K+eij)M  -  (J-K-i+eiJ)M]| 


*  K,je“+  [<j-^vM-  (j-l-i+0-)M 


iJ' 


]h 


(9c) 


K f  L  ~  1 )  2j  •  •  .  ,  p 


and  aM+L,  M+K  aM+K,  M+L' 

Thus,  the  entire  A  matrix  is  symmetric! 

The  elements  of  the  B  matrix  are  given  by 
S* 

b  =  —  2  2  -TT  (J-1+B..f  H=0,  1 . M-l  (9d) 

M.+  1  li!  T  .  Z  lJ 

J  1  0\ 


-ill 

5M+K  M! 


s. 


2  2 
J  i 


6  eM  +  8 

K,  J  iJ  K< J 


[(j-K+eij)M-(j-K-i+eij)M]j 


(9e) 


K  =  1,  2,  . 


where  in  (9a)- (9e)  it  is  understood  that  the  summation  limits  on  J  are  1  to  p 

and  the  limits  on  i  are  N  +1  to  N  . 

J  -  1  J 

Equations  (4),  (8),  and  (9)  are  the  relationships  required  for  ‘'matched11 
polynomial  least-squares  fitting  of  non -polynomial  data  by  means  of  multiple 
M‘ th-degree  polynomials  in  segments,  and  constrained  to  satisfy  continuity 
conditions  through  the  (M-l)  derivative.  The  parameter  estimates  via 
Eqs.  (8)  are  linear  and  the  matrix  elements  are  simple  summations  of  simple 
algebraic  expressions  [Eqs.  (9)  J.  The  special  case  of  ordinary  least -squares 
polynomial  fitting  about  t  =  T^  results  by  setting  p  =  1  in  Eqs.  (5a)  and  (9). 

It  seems  clear  that  computation  time  and  storage  should  be  nearly  minimal 
and  whether  the  technique  has  sufficient  accuracy  compared  to  iterative 
maximum-likelihood  estimation  becomes  the  chief  concern.  Adaptation  to 
real-time  (3  estimation  will  now  be  discussed. 

Appendix  I  particularizes  these  formulas  to  the  case  of  quadratic  poly¬ 
nomials  (M=2).  "Matched  cubics"  (M=3)  are  the  smallest  degree  polynomials 
for  which  the  second  (and  lower  order)  derivatives  are  continuous  throughout 
the  gross  interval  of  interest.  Consequently,  use  of  M=3  in  the  above  formulas 
seems  particularly  promising  for  the  trajectory  applications  that  follow. 
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From  Newtonian  dynamics 


S 


(10) 


where  S  is  the  second  derivative  of  range,  p  the  air  weight  density  at  time 
t  ,  v  the  missile  velocity,  S  the  range  rate,  and  (3  the  ballistic  coefficient. 
Physically,  g/p  is  the  component  of  the  drag  force  along  the  radar  line  of 
sight.  The  f  function  contains  components  of  acceleration  along  the  line  of 
sight  due  to  angular  rates,  gravity,  and  earth  rotation  effects.  No  Doppler 
data  will  be  assumed  in  what  follows.  (Inclusion  of  Doppler  information  will 
be  discussed  later,  ) 

The  f  and  g  functions  in  Eq.  (10)  contain  no  second  derivatives.  The 
smooth  range,  angle  values,  their  rates,  and  S  all  at  time  t  are  to  be  de¬ 
termined  from  the  radar  data  by  means  of  Eqs.  (8)  and  (9).  For  M  =  2,  since 
S  has  different  constant  values  for  each  region,  Eq.  (10)  solved  for  (3  ,  repre 
sents  a  M  stair  case"  ballistic  coefficient  history.  For  M  >  2,  S  is  continuous 
and  consequently,  (3  is  also.  The  following  outlines  a  real-time  P  method. 

1.  A  few  samples,  say  "i^11  ,  of  radar  metric  data  are  ob¬ 
tained  at  a  high  re-entry  altitude.  Times  Tj,  J  =  1,  2,  ... 
have  been  defined  a  priori  and  for  convenience  it  will  be 
assumed  that  at  least  one  radar  sample  will  be  obtained  in 
each  region  (Fig.  1). 

2.  Equations  (8)  and  (9)  with  p  =  1,  for  example,  are  solved 
with  5^  representing  the  range  radar  data  with  c =  ^RANGE^ 

and  i  =  1,  2,  .  .  .  ,  i  .  Thus,  we  obtain  the  first  estimate  of  S .  , 

1  *  1 
...  .  Similarly,  Eqs.  (8)  and  (9),  with  representing 

the  elevation  data  and  then  the  azimuth  data,  are  solved  to  give 

the  first  estimates  of  E^,  E^,  E^  ...  and  then  A^,  A^,  A^... 


^  In  particular,  the  absence  of  angle  second  derivatives  makes  this  equation 
particularly  suitable  for  calculating  the  ballistic  coefficient  because  radar 
angle  errors  are  ordinarily  much  more  serious  than  range  and  Doppler  errors 
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3. 


Equations  (4)  give  S^,  S^,  S^,  E^,  E^,  A ,  A^,  and  Eq.  (10)  can  be 
solved  for  p  (F),  yielding  the  first  estimate  of  |3  in  region  1  at  t  =  t.. 

4.  As  the  next  set  of  one  (or  more)  radar  data  samples  are  re¬ 
ceived,  giving  a  total  of  i^  =  i^  4  Ai^  ,  two  distinct  cases  arise: 

a)  i^  lies  within  the  previous  or  p  =  1  region  above,  or 

b)  i^  lies  in  the  p  ^  2  region. 


For  case  (a)  the  changes  in  the  matrix  elements  Aa^  and  Ab^  ,  a,  (3  = 

1,  2,  .  .  .  ,  (M  4  1)  are  calculated  from  Eqs.  (9)  with  p  =  1  and  the  index 
i  ranging  from  i^  +  1,  i^  .  It  is  important  to  notice  that  only  the  new  data 
samples  are  included  in  these  summations  and  that  the  old  summations  need 
not  be  recalculated.  The  new  matrix  elements,  uv  4  ln  ,  are  given  in  terms 
of  the  old  uv"  as 


a 


(I'+i) 
a  P 


+  Aa 


q(3 


b 


I 

\ 

a 


v  +  1) 


4  Ab 

a 


(11) 


and  the  parameter  estimates  are  updated  via  Eq.  (8a) 

X  =  A-1B  .  (12) 

Equation  (4)  is  then  used  to  calculate  S.,  S.,  S.,  E.,  E.,  A.,  A.,  and  the  new 
v  i  iiiiii 

p  value,  3(F)  ,  is  calculated  from  Eq.  (10). 

For  case  (b),  assume  for  convenience  that  i^  lies  in  the  second  region 
p  =  2  .  The  changes  in  the  matrix  elements  ^aQp  and  AbQ  9  a,  p  =  1,  2,  .  .  .  , 

(M  4  2)  are  calculated  using  only  the  new  radar  data  and  Eqs.  (9).  The  sum¬ 
mation  limits  are  p  =  1  and  i  =  i^  4  1,  for  the  first  region;  and  p  =  2  and 

i  =  4  1,  i^  for  the  second.  The  new  values  of  the  matrix  elements  are 

calculated  from  Eq.  (11)  and  the  updated  parameter  values  from  Eq.  (12). 

Then,  S^,  .  .  .  ,  P(F)  are  calculated  as  for  case  (a)  above.  We  note,  how¬ 

ever,  that  this  case  (b)  not  only  updates  all  of  the  previous  parameter  estimates, 
but  furnishes  the  first  estimates  for  P  in  region  2,  requiring,  however,  the 
inversion  of  a  (M  4  2)  x  (M  4  2)  matrix  instead  of  a  (M  4  1)  x  (M  4  1)  as  in 
case  (a). 
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5,  Beyond  steps  (3)  and  (4),  continuation  of  the  real-time 
calculation  should  be  evident.  For  practicality,  a  limit  on 
the  largest  size  matrix  to  be  inverted  must  be  imposed.  Such 
a  limit  might  in  practice  require  some  additional  steps  which 
need  not  concern  us  here. 

It  was  noted  that  the  A  and  B  matrices  involve  simple  summations 
over  the  available  data  points,  and  as  more  data  points  are  obtained,  these 
summations  need  only  be  updated  rather  than  recalculated.  The  solution 
matrix  X  for  the  latest  parameter  estimates,  however,  do  require  recalcu¬ 
lation.  Attention  is  now  focused  on  this  calculation. 

( y ) 

Assume  that  Xv  1  has  been  determined  and  includes  data  points  recently 
acquired  in  region  J  .  New  data  points  are  then  acquired 

(a)  entirely  in  the  same  region  J  (as  in  case  (a) 
above,  or 

(b)  at  least  one  of  the  new  points  falls  in  region 
(J  +  1)  (as  in  case  (b)  above)  . 

Letting  A^  +  ^,  B^  +  ^,  and  X^  +  ^  denote  the  new  A,  B,  and  X 
matrices  respectively,  we  have 


A("+1> 

b(^d 

x("+1) 


=  +  (AA) 

=  B^)  +  (AB) 
=  X(V)  +  (AX) 


(13) 


The  matrices  with  superscript  (v)  in  Eqs.  (13)  are  known  from  previous 
calculation.  The  (A A)  and  (AB)  matrices  defined  by  Eqs.  (11)  are  updated 
using  Eqs.  (9).  Thus,  only  AX  and  consequently  X^  +  ^  are  to  be  redetermined. 
From  Eqs.  (12)  and  (13)  there  follows 


AX  =  [A^  +  1)]_1[AB  -  (AA)  X*1')]  .  (14) 

Equation  (14)  will  apply  for  both  cases  (a)  and  (b)  providing  X^  ,  for  example, 
is  defined  as  a  ( J  +  1)  +  M  column  vector  with  the  last  element  identically 
zero  for  case  (a)  or  identically  the  (J  +  1  +  M)  independent  variable  if  case  (b); 
similarly  with  the  other  matrices  in  Eq.  (14). 
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The  new  parameter  estimates  are  thus  determined  from  Eq.  (14)  in 

( i> \  li>\ 

terms  of  the  old  parameters  X  ,  the  old  A  matrix  A  ,  and  the  changes 
in  the  A  and  B  matrices.  Computation  time  and  storage  are  clearly  a 
minimum  except  possibly  for  the  inversion  of  the  A^  +  ^  matrix  that  is  re¬ 
quired  with  every  acquisition  of  new  data. 

When  good  Doppler  data  is  available,  the  same  formulas  and  methods 
presented  are  used  with  the  required  changes  in  notation.  The  S(t)  and  S(t) 
histories  previously  calculated  from  the  range  data  are  merely  replaced  by 
the  values  calculated  from  the  Doppler  information.  For  M  =  2  ,  the  (3 
history  calculated  using  Doppler  information  is  clearly  continuous. 

The  A^1  ^  matrix  of  Eq.  (14)  grows  by  one  row  and  one  column  when¬ 
ever  data  is  acquired  in  a  new  region,  which  restricts  the  number  of  regions 
that  can  be  utilized  in  applications  of  this  technique. 

HS:cm 
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APPENDIX  I 


For  convenience,  we  particularize  the  formulas  derived  in  the  text  to 
matched  quadratic  polynomials  and  list  the  results  here  (M  =  2)  . 

Equation  (4)  of  the  text  becomes 

S(ti)  =S1+h(J-l+0  )S1+^9  2S  +h2  E  )Sk  (A-l) 

k=  1 


J-l 


S(ti)  =s1+h(eiJsJ+kZi 


sk> 


e 


iJ 


Tx  *  t.  *  Tt  ,  . 
J  l  J  +  l 


J  =  1  2 

i J  1  “  )  •  •  •  J 


p 


The  matrix  elements  in  conjunction  with  Eqs.  (8a)  and  (8b)  reduce  to 
Nt 

P  J 

a  =  Z  Z  ,  /  2 

11  T  i  -nt  ,  ,  I/O-. 

J  =  1  i=NJ_1  +  l  '  i 

a,  ,=EEh(J  -  1  +  e.T)/o-.2 
J  l 

a22  =  h2  Z  Z  (J  -  1  +  eiJ) 2/cr.2 
J  i 


a2+K,  1  2  9iK  +  6J>K^J'K"2  + 9iJ^/°'i 

J  l 

K  =  1,  2, 

a2+K,2=h3!E(J-1+9iJ>[-TSeiK+6J>K(J-K-7+9iJ>]/Ti 


•  >  P 
2 


J  i 


A-l 


1 2+K,  2+L  h  ^  2  0iK+6J>K(J‘K"  2  +  0iJ^  ^  2  0iL+ 6 J>L(  J _L"  2  +  9iJ) -^i 

J  l 

K  =  1,  2,  .  .  .  ,  p  L  =  1,  2,  ...  p 


aKL  aLK 


Nt 

P  J  >:<  ; 

b ,  =  E  E  S.  /<r. 

1  J=1  i=Nj_  j  +  l  1  1 


b2  =  h  Z  f  (J  -  1  +  eu)  S~/ar.‘ 
J  1 


L,  K  =  1,  2,  .  .  .  ,  2  +  p 


b2  +  K  =  h‘?  ?  eiK  +  6J>K<  J'K4  +  eij' 

J  1 


1 


where 


J>K 


i: 


for  J  >  K 
otherwise 


J,  K 


(  1  for  J  =K 

I  0  otherwise 


(A-2) 


K  =  1,  2,  . 
K 


An  alternative  to  standard  formulas  for  performing  numerical  integration 
is  obtained  by  integrating  Eq.  (A-l).  The  result  is 


I 


:p+i 


2  P 


S(t)dt=hp{Si+S1l|^+s^-  E[3(p-k+l)(p-k)+l]Sk}  ( A-3) 

‘  **  k=  1 

p  —  1 >  ••• 


Referring  to  Fig.  1,  and  S]  are  the  values  of  the  integrand  and  its  first 


1 


derivative  at  t  =  T.  and  S,  is  the  value  of  the  second  derivative  at 

1  k 

k  =  1,  2,  .  .  .  ,  p.  The  formula  becomes  exact  as  p  ->  °°  s  h  -*  0  . 
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Integration  of  Eq.  (5a)  followed  by  summation  over  J  gives  the  more 


general  quadrature  formula 
T 


P+i 

j  S(t)  dt 


M-l 


k+1 


_  V  s  (k)  (M 

Lt  1  (k+1 )! 


M+l 


(M+l ) 


k=0 


k=  1 


sk<M»  X 


|(p-k+  1) 


M+l  .  M+ 1 ) 
-  (p-  k)  f 


p  =  1,  2, 


(A-4) 


The  integrand  S(t)  can  be  evaluated  at  various  points  t^.  The  matrix 
elements  of  Eqs.  (9a-9e)  can  be  calculated,  and  the  independent  parameters 
appearing  in  the  right-hand  side  of  Eq.  (A-4)  are  obtained  via  Eq.  (12). 
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